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EVALUATION  OF  THE  NEAR- WALL  K-6  TURBULENCE 
MODEL  BY  COMPARISON  WITH  DIRECT  SIMULATIONS  OF 
TURBULENT  CHANNEL  FLOW 


I.  INTRODUCTION 

During  the  past  three  decades  the  fluid  dynamics  community  has  focused  much  of 
its  attention  on  the  development  of  techniques  and  capabilities  for  the  calculation  of 
turbulent  flow  fields.  Increases  in  available  computational  capabilities  have  permitted 
researchers  to  develop  and  test  increasingly  sophisticated  and  complete  models  for  the 
numerical  simulation  of  turbulent  flows.  The  state  of  the  art  has  reached  a  point 
where  direct  numerical  simulations  of  turbulent  flows  at  low-to-moderate  Reynolds 
numbers  are  feasible.  Such  work  has  been  reported  in  several  studies  including  those  by 
Kim,  Moin  and  Moser  (1987)  and  Handler,  Hendricks  and  Leighton  (1989),  hereinafter 
referred  to  as  KMM  and  HHL,  respectively.  In  these  simulations,  the  unsteady  Navier- 
Stokes  equations  axe  solved  numerically.  All  essential  scales  of  the  turbulent  flow  are 
resolved  and  no  subgrid  modeling  of  the  turbulence  is  employed.  Such  calculations, 
however,  are  feasible  only  for  research  purposes,  and  calculations  of  turbulent  flows  for 
engineering  purposes  utilize  the  time-averaged  Navier-Stokes  equations  coupled  with 
some  level  of  turbulence  modeling. 

For  many  applications,  the  model  of  choice  is  the  two  equation  k  —  e  model.  This 
model  for  turbulence  transport  came  into  use  in  the  early  1970’s  based  largely  on  the 
work  of  Hanjalic  and  Launder  (1972)  and  Launder,  Reece  and  Rodi  (1975)  and  became 
a  standard  because  of  its  relative  simplicity  and  success  in  providing  good  predictions 
for  a  large  range  of  turbulent  flows.  In  spite  of  some  significant  shortcomings  in  the 
k  —  e  model,  it  still  retains  a  position  as  the  standard  for  comparison. 


Manuscript  approved  May  I.  1989 


1 


In  the  development  of  the  k  —  e  model,  attention  was  focused  on  high  Reynolds 
number  flows,  such  as  thin  shear  layers,  and  the  constants  in  the  model  were  chosen  to 
give  agreement  with  experimental  data.  Consequently,  the  predictive  capabilities  of  the 
model  are  best  for  strongly  turbulent  flows.  For  flows  near  fixed  surfaces,  adjustments 
to  the  model  are  needed.  In  Hanjalic  and  Launder  (1972)  and  Launder,  Reece  and  Rodi 
(1975),  boundary  layer  flows  were  treated  by  matching  the  high  Reynolds  number  flow 
with  wall  functions,  e.g.  “law-of- the- wall”  profiles,  at  the  first  grid  point  away  from  the 
wall.  In  that  way,  the  model  did  not  have  to  be  directly  modified  for  the  low  Reynolds 
number  flow  adjacent  to  the  wall.  In  work  since  then,  Hanjalic  and  Launder  (1976)  and 
a  number  of  researchers  have  proposed  extensions  to  the  basic  high  Reynolds  number 
model  that  would  enable  the  k  —  e  model  to  be  directly  used  all  the  way  to  the  wall. 
Patel,  Rodi  and  Scheuerer  (1985)  review  a  number  of  near-wall  models  and  found  that 
while  the  better  ones  did  quite  well,  there  were  still  areas  where  improvement  was 
needed.  One  of  the  limitations  in  the  development  of  near-wall  models  has  been  the 
difficulty  in  making  accurate  and  reliable  measurements  of  the  flow  quantities  near  the 
wall  and  the  fact  that  some  of  the  correlations  that  are  needed  in  improving  the  k  —  e 
model  cannot  be  obtained  experimentally. 

Data  from  direct  numerical  simulations  of  turbulent  flows  can  be  used  in  place 
of  experimental  data  in  evaluating  and  improving  the  near-wall  turbulence  models. 
Mansour,  Kim  and  Moin  (1988)  have  used  the  data  of  Kim,  Moin  and  Moser  (1987)  to 
analyze  some  of  the  terms  in  the  modeled  transport  equations  for  the  Reynolds  stresses. 
In  this  study,  we  use  the  direct  simulation  data  of  Handler,  Hendricks  and  Leighton 
(1989)  to  examine  the  terms  in  the  k  —  e  model  with  emphasis  on  the  behavior  of  the 
terms  near  the  wall. 

II.  DIRECT  NUMERICAL  SIMULATION 

Numerical  methods  are  being  developed  at  the  U.  S.  Naval  Research  Laboratory 
(NRL)  for  the  direct  numerical  simulation  (DNS)  of  turbulent  flows.  Much  of  the 
background  for  this  effort  is  provided  in  the  study  by  HHL.  In  their  report,  HHL 
review  the  numerical  methods  which  me  currently  in  vogue  and  present  the  results 
of  applying  these  methods  to  the  computation  of  turbulent  channel  flows.  Several 
calculations  parameterized  by  the  extent  of  the  physical  domain,  grid  resolution,  and 
Reynolds  number  are  described.  In  this  study,  we  use  the  Chanl.l  data  of  HHL.  This 
is  one  of  their  lower  resolution  (spatial)  data  sets;  however,  it  is  one  for  which  complete 
velocity  and  pressure  are  available  at  widely  spaced  time  intervals,  thus  insuring  that 
the  separate  realizations  axe  statistically  independent. 
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In  units  of  channel  half-width,  h.  the  channel  dimensions  are  2  x  5  x  5  in  the 
vertical  (12),  the  streamwise  (xi)  and  lateral  (x3)  directions  respectively.  The  flow  is 
computed  with  16  x  64  evenly  spaced  grid  points  in  the  horizontal  (ij  —  x3)  plane  and 
33  points  with  Chebyshev  scaling  are  used  in  the  X2  direction.  The  Reynolds  number, 
Re,  is  2215  based  on  h  and  Uq,  the  initial  laminar  centerline  velocity.  The  governing 
equations  axe  nondimensionalized  by  the  wall  shear  velocity,  a  viscous  length  unit  and 
a  viscous  time  unit.  The  wall  shear  velocity,  ur,  is  given  by 
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The  viscous  length  unit  (or  wall  unit)  is  lv  =  v/uT  and  the  viscous  time  unit  is  tv  = 
v/u In  the  nondimensionalized  units,  the  Reynolds  number  based  on  the  wall  shear 
velocity  is  125.  This  is  lower  than  the  wall  Reynolds  number  of  ISO  for  the  KMM  data 
and  the  grid  resolution  of  the  Cham  1.1  data  is  lower.  HHL  conclude,  however,  that 
the  mean  statistics  of  their  data  set  compare  well  with  experiment  and  with  the  more 
highly  resolved  data  of  KMM  insofar  as  one-point  statistics  are  concerned. 

The  equations  are  solved  using  a  pseudospectral  method  in  which  Chebyshev  ex¬ 
pansion  is  used  in  the  wall  normal  (x2)  direction  and  Fourier  series  are  used  in  the 
streamwise  (xi)  and  spanwise  (x3)  directions.  The  domain  size  in  wall  units  is  640, 
250  and  640  in  the  Xj,  x2  and  x3  directions,  respectively.  In  the  x2  direction,  the  use 
of  Chebyshev  polynomials  gives  a  variable  step  size,  with  Ax2  «  0.6  at  the  wall  and 
Ax2  «  12  at  the  channel  center.  In  the  Xj  and  x3  directions  the  step  size  is  uniform 
and  is  Axj  «  40  and  Ax3  ~  10.  In  the  direct  numerical  simulation  used  in  the  present 
work,  there  are  33  distinct  realizations  of  the  velocity  and  pressure  fields  which  are 
separated  by  50  viscous  time  units.  This  separation  in  time  is  sufficient  to  insure  that 
each  realization  of  the  fields  is  statistically  independent. 

We  denote  the  instantaneous  values  of  variables  by  uppercase  symbols,  the  time 
averaged  values  with  an  overbar,  and  the  fluctuations  from  the  time  averages  by  low¬ 
ercase  symbols  so  that,  for  example, 

Ui(t)  =  Ui  +  u,(t)  and  P(t)  =  P  +  p(t). 
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Also,  we  denote  the  turbulence  intensities  as  u'  =  y  u-.  The  velocity  field  in  the 
channel  is  periodic  in  the  £3  and  £3  directions  so  that  averaging  over  the  x\  —  £3  plane 
at  each  value  of  x%  gives  dependent  variable  profiles  which  are  functions  of  the  wall 
normal  coordinate  (£2).  Averages  are  obtained  by  summing  the  individual  profiles 
over  33  separate  realizations.  Correlations  of  the  fluctuating  variables  such  as  u\U2  are 
computed  by  forming  the  product  at  each  grid  point  in  the  domain,  averaging  over  the 
£1  —  £3  plane,  and  then  averaging  over  the  33  realizations. 

The  primary  averaged  data  from  the  direct  numerical  simulation  of  channel  flow 
by  HHL  are  shown  in  Figs.  1  and  2.  Figure  1  shows  the  mean  streamwise  velocity  U\ 
and  the  root  mean  square  of  the  pressure  fluctuations  p' .  The  velocity  profile  shows  a 
fully  developed  turbulent  flow  profile  with  channel  center  value  of  17.  The  rms  of  the 
pressure  fluctuations  show  peaks  of  1.9  at  y/h  =  ±100.  At  the  channel  center,  p'  has  a 
minimum  value  less  that  1.0.  The  values  at  the  walls  are  between  1.2  and  1.3.  The  wall 
values  are  lower  than  those  reported  by  KMM  (1.5),  and  the  difference  is  attributable 
to  the  lower  Reynolds  number  of  the  HHL  simulation.  Figure  2  shows  the  turbulence 
intensities,  u\,  u'2  and  u'3 ,  and  the  Reynolds  shear  stress,  u  1  ±2 .  It  is  noted  by  HHL 
that  the  location  of  the  peak  values  of  u[  are  in  excellent  agreement  with  the  data  of 
KMM  and  experimental  data,  although  the  peak  values  are  slightly  lower.  The  overall 
agreement  between  the  HHL  data  and  KMM  data  is  good.  The  Reynolds  shear  stress, 
«itt2,  shows  peak  values  of  ±0.65  at  y/h  =  ±95  and  linear  variation  over  the  interval 
—85  <  y/h  <  85.  The  averaged  data  in  Figs.  1  and  2  exhibit  good  overall  symmetry 
between  the  channel  walls. 

III.  THE  k  -  e  TRANSPORT  EQUATIONS 

Exact  Transport  Equations 

The  two  equation  k  —  e  model  for  turbulent  flow  requires  transport  equations  for 
the  turbulence  kinetic  energy,  k ,  and  for  e,  the  rate  of  dissipation  of  turbulence  energy. 
The  exact  equation  for  the  transport  of  turbulence  kinetic  energy  is  obtained  from  the 
sum  of  the  transport  equations  for  the  three  normal  Reynolds  stresses.  With  a  slight 
rearrangement,  the  k  transport  equation  in  indicial  notation  as  given  by  Hinze  (1975) 
is, 
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(1) 


(2) 


(3) 


—  v 

(4)  (5) 

where  k'  is  half  the  sum  of  the  normal  stresses  and  k  is  its  averaged  value.  It  should 
be  noted  here  that  for  the  channel  flow  simulations  the  mean  flow  is  homogeneous  in 
the  horizontal  coordinate  directions  so  that  U\  —  f(x 2)  and  U2  =  U3  =  0.  Because 
of  symmetry  with  respect  to  planes  normal  to  the  spanwise  direction,  all  correlations 
involving  u3  and  uneven  derivatives  with  respect  to  x3  are  zero. 


dui  duj\  duj 
dxj  dx i  )  dx{ 


(1) 


d  ( dui  duj 
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The  terms  on  the  right  hand  side  correspond  to  (1)  production  of  turbulence  kinetic 
energy,  (2)  transport  of  turbulence  energy  (3)  transport  of  the  pressure  fluctuations, 
(4)  viscous  diffusion  and  (5)  viscous  dissipation.  For  an  incompressible  fluid,  terms  (4) 
and  (5)  can  be  written  in  the  alternate  form, 


Dk 

Dt 
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(5) 
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While  the  sum  of  terms  (4)  and  (5)  is  the  same  in  Eqs.  1  and  2,  the  individual 
terms  are  not  equivalent.  Term  (5)  in  Eq.  1  is  the  full  dissipation  term  and  term  (5) 
in  Eq.  2  is  the  dissipation  only  in  homogeneous  turbulence.  It  is  the  latter  term  which 
is  modeled  in  the  k  —  e  transport  equations. 

The  exact  transport  equation  for  the  dissipation  rate  term  is  derived  from  the 
unsteady  Navier-Stokes  equations.  We  use  the  form  given  by  Hanjalic  and  Launder 
(1976), 


(1)  (2) 
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The  first  term  on  the  right  hand  side  of  Eq.  3  represents  turbulent  production  of  e. 
Term  (2)  represents  dissipation  of  c.  The  terms  in  the  first  pair  of  braces,  terms  (3),  (4) 
and  (5),  represent  turbulent  transport,  pressure  transport  and  viscous  diffusion.  The 
last  three  terms,  terms  (6),  (7)  and  (8),  are  production  terms. 

The  individual  terms  in  Eqs.  1-3  have  been  calculated  from  the  DNS  data  of  HHL. 
Products  such  as  u,Uj  are  formed  at  each  grid  point  in  the  domain,  the  correlation  is 
then  averaged  over  the  xi  -  X3  plane  giving  a  profile  in  X2.  The  X2  profiles  are  then 
averaged  over  the  33  realizations  in  the  saved  data  sets.  In  forming  the  averages  of  the 
terms,  we  make  use  of  the  symmetry  (in  the  mean)  of  the  channel  flow  about  the  center 
line.  In  the  following  figures,  we  use  the  wall  normal  distance  y+  with  origin  at  the 
wall  whereas  in  the  preceding  figures  we  used  y/h  with  origin  at  the  channel  center. 

Figure  3  shows  the  sums  of  terms  (4)  and  (5)  in  Eqs.  1  and  2.  The  agreement 
between  the  sums  of  these  terms  in  nearly  exact,  an  indication  of  the  fidelity  of  the 
simulation  in  maintaining  a  divergence  free  flow.  Comparison  of  terms  (5)  in  the  two 
equations  shows  only  slightly  greater  differences  than  is  seen  in  the  sums  of  the  terms. 
The  individual  terms  in  Eq.  2  are  shown  in  Fig.  4.  These  data  agree  very  well  with  the 
data  presented  by  Mansour,  Kim  and  Moin  (1988),  hereinafter  referred  to  a s  MKM, 
with  the  differences  (small)  being  attributable  primarily  to  the  differences  in  resolution 
of  the  two  calculations  and  to  a  lesser  extent  to  the  differences  in  Reynolds  numbers. 
Both  HHL  and  KMM  make  extensive  comparisons  with  experimental  data  and  find 
very  good  agreement,  and  such  comparisons  are  not  included  herein.  The  interested 
reader  may  refer  to  those  sources. 

As  seen  from  Fig.  5,  across  most  of  the  channel,  the  principal  contribution  to  the 
transport  of  e  comes  from  terms  (1)  and  (2).  Near  the  wall  (for  y+  <  25)  terms  (6) 
and  (7)  become  significant  contributors  and  very  near  the  wall  (for  y+  <  3),  term  (5) 
becomes  a  major  contributor.  At  the  wall,  it  is  term  ^5)  which  balances  term  (2). 
The  other  three  terms  ((3),  (4)  and  (8))  contribute  very  little  to  the  transport  of  e  in 
the  channel  and  are  not  shown  in  this  figure.  Except  for  the  small  differences  due  to 
resolution  and  Reynolds  number  differences,  these  data  agree  nearly  exactly  with  the 
corresponding  data  of  MKM. 

For  engineering  calculations  of  turbulent  flows,  the  k  —  e  model  is  adequate  for 
many  applications,  but  the  terms  involving  the  fluctuating  velocity  and  pressure  com¬ 
ponents  must  be  modeled.  The  models  which  are  in  use  have  been  developed  based  on 
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order  of  magnitude  analysis  of  the  transport  equations  and  contain  empirical  constants 
which  have  been  determined  from  experimental  data.  The  experimental  data  base  is 
incomplete  because  of  difficulties  of  measuring  some  of  the  terms  in  the  exact  equations. 
The  data  from  the  direct  simulations  permit  a  precise  examination  of  the  terms  in  the 
transport  equations  and  a  new  evaluation  of  the  empirical  constants  in  the  turbulence 
models. 

Modeled  Transport  Equations 


In  modeling  the  terms  of  the  k  transport  equations,  terms  (2)  and  (3)  are  taken 
together  as  the  diffusional  flux  of  k  and  and  assumed  proportional  to  the  gradient  of 
k.  Rodi  (1980)  gives 

**,■(*  +  £)  (4) 

\  PJ  Ok  OX  I 

with  ak  an  empirical  diffusion  constant.  Term  (4)  of  Eq.  2  is  easily  combined  with 
Eq.  4,  and  term  (5)  is  simply  e.  The  Reynolds  stresses  u,u}  are  modeled  using  the 
Boussinesq  eddy  viscosity  relation 
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with  6tJ  the  Kronecker  delta.  The  eddy  viscosity,  turbulence  kinetic  energy  and  the 
rate  of  dissipation  of  turbulence  kinetic  energy  are  related  by 

jfc 2 

(6) 

where  Cy  is  a  constant.  This  expression  for  vt  yields  good  agreement  for  flows  such 
as  thin  shear  layers,  but  for  turbulent  flow  next  to  a  wall,  C y  is  not  constant  or  a 
wall  damping  function  needs  to  be  included.  Patel,  Rodi  and  Scheuerer  (1985)  use  the 
model 

vt-CyfJ^j,  (7) 

with  the  functional  form  of  fy  to  be  determined  by  comparison  with  experimental  data. 
The  modeled  form  of  the  k  transport  equation  then  becomes 
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The  steps  in  developing  the  model  for  the  e  transport  equations  are  quite  complex 
and  are  not  detailed  here.  Hanjalic  and  Launder  (1976)  give  the  following  form  which 
retains  a  direct  Reynolds  stress  expression 
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They  note  that  at  a  wall,  the  e2 /k  term  is  unbounded,  since  e  has  a  finite  value  at  the 
wall  and  k  — ♦  0.  They  eliminate  this  difficulty  by  replacing  e2  by  ee  where  e  is  defined 

by 
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and  shew  that  e/k  approaches  a  finite  value  at  the  wall  which  ensures  that  ee/k  is 
bounded  as  y+  — ♦  0. 

For  the  k  —  e  model  equations,  the  last  term  in  Eq.  10  is  written  as  a  diffusion 
term  with  UiUjk/e  replaced  by  —  iq/cre,  where  a,  is  an  empirical  dissipation  constant. 
This  term  is  then  combined  with  term  (5)  of  Eq.  3.  The  modeled  transport  equation 
for  e  takes  the  form 
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with  ut  and  V  as  given  above.  Commonly  accepted  values  of  the  constants  ak,  ot,  C( \ 
and  C( 2  are  1.0,  1.3,  1.44  and  1.92,  respectively. 

For  the  channel,  U 2  and  U 3  are  zero  as  are  the  derivatives  of  the  averaged  tur¬ 
bulence  variables  in  the  rj  and  13  directions.  The  modeled  transport  equations  thus 
reduce  to  the  following  form 


Dk 

~Dt 


d 

dx  2 


dk 
dx  2 


+  V-e 


(13) 


and 


8 


Dt 

Dt 


d 

dx2 


dt  ' 
dx2  _ 


+c-'iv~c4 


(14) 


with 
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The  eddy  viscosity  is  given  by  Eq.  7.  For  thin  shear  flows  C ^  is  a  constant  and  and 
a  value  of  0.09  is  typical.  Likewise,  for  thin  shear  flows,  f ^  =  1.0.  For  wall  bounded 
flows,  is  the  wall  damping  function  and  is  usually  chosen  to  approach  1.0  when  the 
flow  is  not  influenced  by  the  wall. 


IV.  EVALUATION  OF  THE  k  -  e  MODEL 


One  of  the  significant  differences  between  the  k  —  t  turbulence  model  and  models 
which  require  additional  transport  equations  is  in  the  treatment  on  the  Reynolds  stress 
terms,  ujuj.  In  the  models  using  separate  equations  for  the  transport  of  the  Reynolds 
stresses,  the  u,ity  terms  are  solved  as  unknown  variables.  In  the  k  —  t  model,  the 
Reynolds  stress  terms  are  modeled  using  the  eddy  viscosity  as  given  in  Eqs.  7  and  16. 
The  accuracy  of  this  model  is  then  central  in  considering  the  accuracy  of  the  k  —  e 
model 

If  ffi  is  taken  to  be  1.0,  the  k  —  e  model  corresponds  to  a  high  Reynolds  number 
model.  Using  the  HHL  data,  we  calculate  the  terms  in  each  side  of  Eq.  16  with  ut  given 
by  Eq.  7,  =  0.09  and  =  1.0.  Figure  6  displays  the  two  sides  of  Eq.  16  and  we  see 

that  the  equation  applies  very  well  in  the  center  of  the  channel  but  that  in  the  near¬ 
wall  region  the  modeled  term  is  much  larger  than  the  Reynolds  stress  computed  from 
the  DNS  data.  "  Vse  det^  show  that  a  damping  term  is  required  to  bring  the  modeled 
Reynolds  str°  u  »to  better  agreement  with  that  calculated  from  the  simulation  data. 

A  number  ot  wall  damping  models  have  been  proposed.  One  of  the  simplest  (and 
oldest)  is  that  p:  ,  ed  by  van  Driest  (1954)  in  which  the  damping  depends  only  on 
the  law-of-  hc-v.vul  coordinate,  y+,  and  a  constant  of  the  turbulence,  A+ .  In  that 
model,  the  appropriate  damping  function  for  the  eddy  viscosity  is  determined  to  be 
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(17) 


fv  ~  1 -exp(-y+/A+)  . 

This  model  has  given  excellent  results  for  simple  flows  like  flat  plates,  but  has  been 
judged  inadequate  for  more  complex  flows,  such  as  those  with  recirculation  and  those 
with  pressure  gradients.  Many  of  the  proposed  alternatives  to  Eq.  17  depend  upon  a 
Reynolds  number  of  turbulence,  either  Ry  =  '/ky+ ,  or  Rt  =  k2 /e,  and  do  not  depend 
solely  upon  y+. 

In  the  review  of  Patel,  Rodi  and  Scheuerer  (1985),  eight  different  models  for  cal¬ 
culation  of  near-wall  turbulent  flows  are  discussed.  They  found  that  the  model  of  Lam 
and  Bremhorst  (1981)  provided  the  most  reasonable  wall  damping  for  the  experimental 
data  they  considered.  In  that  model,  the  damping  function  is  given  by, 

U  =  1  —  exp(— 0.0165  Ry)  +  (18) 

Patel  and  his  colleagues  have  since  used  the  Lam-Bremhorst  model  for  a  number  of 
turbulent  flow  calculations  with  good  results  (e.g.  Chen  and  Patel,  1987). 

In  the  review  of  Patel,  et  al.,  C ^  =  0.09  was  used  for  most  of  the  models,  including 
the  Lam-Bremhorst  model.  Since  0.09  is  a  commonly  accepted  value  for  C„,  we  use 
that  value  as  the  basis  for  comparison.  Figure  7  shows  the  Reynolds  shear  stress  as 
calculated  from  the  DNS  data  and  as  modeled  by  Eqs.  16  and  7  using  CM  =  0.09  and 
f n  from  Eq.  18.  With  this  wall  damping  the  modeled  and  exact  shears  have  about  the 
same  peak  value,  but  at  different  distances  from  the  wall.  However,  the  model  gives 
too  low  a  value  for  y+  >25. 

Figures  8  and  9  show  respectively  the  budgets  for  the  terms  in  the  modeled  k  and 
e  transport  equations  using  the  wall  damping  function  in  Eq.  18.  For  the  k  transport 
equation,  the  diffusion  term  is  very  small  in  the  center  of  the  channel,  then  becomes 
more  significant  for  y+  <  30,  peaks  at  y+  =  15,  changes  sign  and  is  large  at  the 
wall  (due  to  the  molecular  viscosity).  The  production  term  is  small  in  the  center  of 
the  channel,  starts  to  increase  for  y+  <  50,  peaks  at  a  large  value  at  y+  =  15  and 
goes  to  zero  at  the  wall.  The  dissipation  term  is  larger  (with  opposite  sign)  than  the 
production  term  for  50  <  y+  <  125,  increases  in  magnitude  as  the  wall  is  approached, 
and  reaches  its  minimum  at  the  wall,  where  it  is  balanced  by  the  diffusion  term.  Also 
shown  is  the  sum  of  the  three  modeled  terms,  and  it  is  seen  that  there  is  significant 
imbalance.  The  source  of  imbalance  is  the  contribution  of  the  production  term.  For 
y+  <  25,  the  production  term  contributes  too  much,  and,  for  y+  >  25,  the  production 
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term  contributes  too  little.  For  the  e  transport  equation,  the  character  of  the  terms  is 
like  that  of  the  terms  of  the  k  transport  equation. 

While  good  results  have  been  obtained  with  the  Lam-Bremhorst  wall  damping 
model  (Chen  and  Patel,  1987),  the  amount  of  imbalance  in  the  modeled  forms  of  the 
transport  equations  suggests  that  the  production  term  can  be  better  modeled.  We  use 
the  HHL  data  to  determine  the  form  that  should  take  and  compare  that  with  the 
form  given  by  the  Lam-Bremhorst  model  and  other  models.  Equations  7  and  16  are 
combined  to  give 
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with  to  be  chosen  so  that  goes  to  1.0  at  the  channel  center.  We  note  that  at  the 
wall  k  =  0  and  at  the  channel  center  dU  \/dx2  =  0.  Therefore,  we  limit  evaluation  of 
Eq.  19  to  the  interval  0  <  y+  <  125. 

Figure  10  shows  the  distribution  of  /M,  across  the  half-channel,  calculated  from 
Eq.  19  and  with  C M  =  0.115,  a  value  that  allows  /M  to  approach  1.0  at  the  channel 
center.  Also  shown  is  /M  from  Eq.  18.  For  10  <  y+  <  30,  the  model  gives  too  large  a 
value  by  up  to  about  70%.  For  y+  >  30,  the  model  gives  as  much  as  one  fourth  too 
low  a  value.  We  have  examined  how  the  values  of  the  model  change  as  the  coefficients 
of  Rt  and  Ry  are  changed.  Since  the  second  term  of  Eq.  18  is  greater  than  one  for  all 
values  of  Rt ,  increases  in  the  coefficient  of  Ry  can  give  values  of  greater  than  one 
at  the  channel  center.  We  find  that  the  coefficients  given  in  Eq.  18  seem  to  give  the 
best  overall  agreement  between  the  model  and  f ^  from  the  DNS  data.  What  is  needed 
for  better  agreement  with  the  HHL  data  is  a  model  which  gives  larger  values  of  for 
the  middle  portion  of  the  channel  and  lower  values  near  the  wall.  We  have  looked  at  a 
number  of  other  wall  damping  models,  including  those  considered  by  Patel,  Rodi  and 
Scheuerer  (1985),  that  of  Jones  and  Launder  (1972)  and  that  of  Kim  (1988).  While 
we  do  not  find  good  agreement  between  the  models  based  on  Rt  and/or  Ry  and  /„ 
calculated  from  the  HHL  data,  the  Lam-Bremhorst  model  does  give  better  agreement 
with  the  HHL  data  than  do  the  other  models. 

Figure  10  shows  a  characteristic  of  which  is  seldom  noted.  As  calculated  from 
the  DNS  data  using  Eq.  19,  /M  has  a  minimum  value  of  0.04  at  y+  =  6  and  then 
increases  toward  the  wall.  Such  a  characteristic  also  appears  in  the  /M  which  Patel, 
et  al.  derived  from  experimental  data  and  show  in  their  Fig.  2.  However,  they  do  not 
comment  on  this  characteristic.  We  have  reproduced  as  the  full  rectangles  in  Fig.  10 
some  of  the  data  from  their  Fig.  2.  The  experimental  data  for  are  in  good  agreement 


11 


with  the  HHL  data  and  show  a  non-zero  minimum  at  about  y+  =  7.  Expansion  of  the 
terms  of  Eq.  19  in  Taylor  series  about  y+  =  0  shows  that  for  small  y+,  itj  =  ay+  +  . . .  , 
U1U2  — >  6y+3  -f-  . . .  ,  e'  —>  c  +  dy+  +  .  • .  and  k'  — *  ey+2  +  . . .  .  Since  near  the  wall 
dU  1/6x2  is  constant,  the  damping  function  in  Eq.  19  must  vary  as  y/y+  +  h  +  .  .  . 
The  result  of  this  simple  analysis  is  difficult  to  work  with  since  it  contains  a  singularity, 
but  it  does  indicate  that  /M  should  have  a  minimum  located  off  the  wall.  None  of  the 
wall  damping  models  vary  as  l/y+  for  small  y+  and  most  go  to  zero  instead  of  to  a 
non-zero  minimum. 

During  the  course  of  examining  the  wall  damping  models  mentioned  above,  we 
also  compared  the  model  of  van  Driest  (1954)  given  in  Eq.  17  with  the  DNS  data, 
and  we  find  that  the  van  Driest  model  agrees  better  with  the  HHL  data  than  does 
the  Lam-Bremhorst  model.  Figure  11  shows  /M  from  the  HHL  data  and  from  Eq.  17 
with  A+  =  26  which  is  the  commonly  accepted  value.  We  note  that  the  model  gives 
too  large  a  value  in  the  interval  10  <  y+  <  50  but  that  the  shape  of  the  distribution 
of  fM  from  the  model  and  from  the  HHL  data  is  similar.  By  introducing  an  effective 
origin  for  y+,  the  modeled  and  that  from  the  HHL  data  can  be  brought  into  very 
close  agreement.  The  use  of  an  effective  origin  mimics  the  observed  behavior  that 
f n  possesses  a  minimum  off  the  wall,  however  it  does  not  produce  a  l/y+  variation. 
Figure  12  shows  a  comparison  between  fu  from  the  HHL  data  and  from  Eq.  17  with 
y+  replaced  by  y+  —  yf~  and  y£  =  6.  With  this  addition  of  an  effective  origin,  the 
agreement  between  from  the  data  and  from  the  van  Driest  model  is  good  over  all 
and  is  almost  exact  in  the  neighborhood  of  y+  =  25.  Near  the  wall,  the  model  gives  a 
value  of  zero  at  the  effective  origin,  whereas  the  HHL  data  show  a  value  of  about  0.04. 
The  model  can  be  brought  into  closer  agreement  with  the  HHL  data  by  introducing  a 
minimum  value,  fo.  The  modified  damping  function  is  then  given  by 

/M  =  /o  +  (1  -  /o)(l  -  exp[-(y+  -  y^)M+])  •  (20) 

The  introduction  of  fo  shifts  the  curve  upwards  and  yj"  needs  to  be  correspondingly 
increased  to  maintain  agreement  with  the  direct  simulation  data.  Figure  13  shows  a 
comparison  between  /M  from  the  HHL  data  and  from  Eq.  20  with  y$=  8  and  f0  =  0.04. 
The  agreement  is  seen  to  be  very  good.  Also  shown  in  this  figure  are  the  Patel,  et  al. 
data  presented  in  Fig.  10.  The  modifications  to  the  van  Driest  wall  damping  function 
are  consistent  with  the  experimental  data  and  the  direct  simulation  data. 

With  f n  as  shown  in  Fig.  13,  and  with  CM  =  0.115,  the  modeled  shear  agrees  very 
well  with  the  exact  shear  as  is  shown  in  Fig.  14.  Also  shown  in  Fig.  14  is  the  shear 
calculated  with  the  unmodified  van  Driest  damping  function.  With  the  unmodified 
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damping  function,  CM  is  taken  to  be  0.09  for  direct  comparison  with  the  modeled  shear 
calculated  with  the  Lam-Bremhorst  model  shown  in  Fig.  7.  That  comparison  shows 
that  the  van  Driest  model  gives  a  greater  shear  across  all  the  channel.  For  y+  <  25, 
the  shear  modeled  with  the  Lam-Bremhorst  damping  function  agrees  better  with  the 
exact  shear  than  does  that  with  the  Driest  damping  function.  However,  for  y+  >  25, 
the  van  Driest  damping  function  yields  the  better  agreement.  Overall,  the  van  Driest 
model  gives  better  agreement  with  the  exact  shear  from  the  HHL  data  than  does  the 
Lam-Bremhorst  model. 

Use  of  the  modified  van  Driest  model  for  reduces  the  imbalance  in  the  sum 
of  the  modeled  terms  in  the  k  transport  equation.  Figure  15  shows  the  variation  of 
the  terms  of  the  k  transport  equation  modeled  with  the  modified  van  Driest  damping 
function.  Comparison  with  Fig.  8  shows  that  the  imbalance  is  reduced  by  more  than  a 
factor  of  two. 

The  DNS  data  also  make  it  possible  to  evaluate  the  validity  of  one  of  the  basic 
assumptions  of  developing  the  k  —  t  model:  that  the  production  and  dissipation  of 
turbulence  kinetic  energy  are  in  approximate  balance.  The  production  of  turbulence 
kinetic  energy  is  calculated  from  the  DNS  data  using  the  equation 

P  =  (21) 

The  variation  across  the  channel  of  the  ratio  V/t  is  shown  in  Fig.  16.  Except  for 
10  <  y+  <  35,  the  ratio  is  less  than  one  and  is  zero  at  the  wall  and  the  channel  center. 
At  y+  =  15,  V/t  has  its  peak  value  of  approximately  1.9.  Averaged  across  the  channel, 
V  / 1  is  less  than  one. 

Rodi  (1980)  reports  the  results  of  considering  a  number  of  turbulent  flows  for 
which  the  ratio  V /t  was  not  unity.  In  earlier  work,  he  had  found  that  the  experimental 
data  correlated  with  a  function  C ^  =  f(V/t)  whereby  C ^  varied  inversely  with  V/t. 
Denoting  by  C^y)  the  product  of  and  /M,  we  show  in  Fig.  17  the  values  of  C^(y) 
plotted  vs.  V/t  with  the  omission  of  the  5  points  nearest  the  wall.  This  shows  a 
distribution  of  CM(y)  vs  V/t  much  like  that  given  by  Rodi  except  that  C^y)  is  constant 
at  about  0.11  for  V/t  <  0.7  instead  of  continuing  to  increase  as  V/t  —*  0. 

V.  SUMMARY  AND  CONCLUSIONS 

In  this  report,  we  have  performed  an  evaluation  of  the  k  —  t  turbulence  model  using 
the  direct  numerical  simulation  data  of  Handler,  Hendricks  and  Leighton  for  turbulent 
channel  flow.  With  the  DNS  data,  we  have  evaluated  the  terms  of  the  exact  k  —  t 


transport  equations  and  have  found  nearly  exact  agreement  with  the  budget  data  pre¬ 
sented  by  Mansour,  Kim  and  Moin.  Evaluation  of  the  terms  of  the  modeled  transport 
equations  shows  the  necessity  of  using  a  wall  damping  function  in  the  modeling  of  the 
eddy  viscosity.  From  the  HHL  direct  simulation  data,  we  have  determined  the  shape 
required  for  a  damping  function  to  obtain  agreement  between  the  modeled  Reynolds 
shear  stress  and  that  calculated  from  the  DNS  data.  We  have  shown  that  the  well  ac¬ 
cepted  Lam-Bremhorst  wall  damping  model  gives  reasonable  agreement  with  the  DNS 
data. 

We  have  also  examined  a  number  of  other  wall  damping  models.  Those  which 
are  functions  of  the  turbulence  Reynolds  numbers  generally  agree  less  well  with  the 
DNS  function  than  does  the  Lam-Bremhorst  model.  We  find  that  the  standard  van 
Driest  wall  damping  model  gives  better  agreement  with  the  DNS  data  than  the  Lam- 
Bremhorst  model.  We  further  find  that  modifying  the  van  Driest  damping  model  to 
include  an  effective  origin  gives  very  good  agreement  between  the  modeled  shear  and 
the  shear  calculated  from  the  HHL  data. 
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Fig.  1  —  Profiles  of  mean  streamwise  velocity,  U) ,  and  root  mean  square  of  pressure  fluctuations,  p\  from  the 
direct  simulation  data  of  Handler,  Hendricks  and  Leighton  (i989). 
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Fig.  3  —  Comparison  of  the  sum  of  the  viscous  terms  in  the  two  forms  of  the  exact  k  transport  equation. 
Terms  (4)  and  (5)  of  Eq.  1  are  shown  as  (+).  Terms  (4)  and  (5)  of  Eq.  2  are  shown  as  (A). 
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Fig.  5  —  Budgets  of  the  principal  terms  of  the  exact  e  transport  equation.  Turbulent  production  of  t  (O);  dis¬ 
sipation  of  t  (  +  );  viscous  diffusion  (*);  production  terms;  term  (6),  (O),  and  term  (7),  ( x). 
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Fig.  12  —  Comparison  of  the  wall  damping  functi 
C„  =  0.115,  (O),  and  from  the  modified  van  Driest  it 
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Fig.  13  —  Comparison  of  the  wall  damping  function,  from  the  experimental  data  of  Patel,  et  al.  (□),  calcu¬ 
lated  from  the  HHL  data,  Eq.  19  with  =  0.115,  (O),  and  from  the  modified  van  Driest  model,  Eq.  20, 
with  >o  =  8  and/0  =  0.04,  (  +  ). 
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Fig.  14  —  Comparison  of  Reynolds  shear  stress,  U!,u2,  calculated  from  the  HHL  data  (O),  modeled  using  the 
modified  van  Driest  wall  damping  model,  Eq.  20,  with  f0  =  0.04,  yo  =  8  and  C„  =  0.115,  (*),  and  modeled 
with  the  standard  van  Driest  wall  damping  model,  Eq.  17,  with  =  0.09,  (+). 
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Fig.  15  —  Budgets  of  the  terms  in  the  modeled  transport  equation  for  turbulence  kinetic  energy  using  the  modi 
fied  van  Driest  wall  damping  model,  Eq.  20,  with  f0  =  0.04,  y0+  =  8  and  C„  =  0.115.  Dissipation  term,  e 
(A);  diffusion  term  (first  term  in  Eq.  13)  (+);  production  term,  P,  (O);  and  sum  of  terms  (□). 
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Fig.  17  —  The  eddy  viscosity  coefficient,  C(y),  as  a  function  of  the  ratio  of  the  production  of  turbulence  kinetic 
energy,  P,  to  the  dissipation  rate,  e. 
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